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ABSTRACT 

A  two-step  hybrid  perturbation-Galerkin  technique  for  improving  the  usefulness  of  per¬ 
turbation  solutions  to  partial  differential  equations  which  contain  a  parameter  is  presented 
and  discussed.  In  the  first  step  of  the  method,  the  leading  terms  in  the  asymptotic  expan¬ 
sion^)  of  the  solution  about  one  or  more  values  of  the  perturbation  parameter  are  obtained 
using  standard  perturbation  methods.  In  the  second  step,  the  perturbation  functions  ob¬ 
tained  in  the  first  step  are  used  as  trial  functions  in  a  Bubnov-Galerkin  approximation.  This 
semi-analytical,  semi-numerical  hybrid  technique  appears  to  overcome  some  of  the  draw¬ 
backs  of  the  perturbation  and  Galerkin  methods  when  they  are  applied  by  themselves,  while 
combining  some  of  the  good  features  of  each.  The  technique  is  illustrated  first  by  a  simple 
example.  It  is  then  applied  to  the  problem  of  determining  the  flow  of  a  slightly  compressible 
fluid  past  a  circular  cylinder  and  to  the  problem  of  determining  the  shape  of  a  free  surface 
due  to  a  sink  above  the  surface.  Solutions  obtained  by  the  hybrid  method  are  compared  with 
other  approximate  solutions,  and  its  possible  application  to  certain  problems  associated  with 
domain  decomposition  is  discussed. 
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1.  Introduction 


A  two-step  hybrid  analysis  technique,  which  combines  perturbation  techniques  with  the 
Galerkin  method,  has  been  presented  and  discussed  by  the  authors.  It  was  applied  to  some 
singular  perturbation  problems  in  slender  body  theory  [5],  as  well  as  to  several  classes  of 
problems  involving  ordinary  differential  equations  [2,6,7].  That  technique  also  promises  to 
be  useful  in  the  analysis  of  a  very  wide  variety  of  partial  differential  equation  type  problems 
as  well.  In  this  paper  we  apply  the  method  to  some  problems  involving  several  independent 
variables.  In  particular,  we  demonstrate  its  usefulness  by  applying  it  to  a  linear  boundary 
value  problem,  a  nonlinear  boundary  value  problem,  and  to  a  (nonlinear)  free  boundary 
value  problem. 

The  method  is  based  upon  a  hybrid  technique  which  was  apparently  first  studied  by 
Ahmed  K.  Noor  and  collaborators  in  conjunction  with  the  finite  element  analysis  of  geomet¬ 
rically  nonlinear  problems  in  structural  mechanics  (see  Geer  and  Andersen  [5]  for  several 
references).  The  Galerkin  method  has,  of  course,  been  known  and  used  for  a  long  time. 
However,  a  principle  problem  associated  with  its  successful  application  lies  in  the  choice  of 
appropriate  basis  functions.  In  a  series  of  papers  Noor  and  his  collaborators  have  shown  for 
a  variety  of  structural  mechanics  problems  that  the  first  few  terms  in  a  Taylor  series  expan¬ 
sion  of  the  solution  of  a  parameterized  system  of  discretized  equations  can  be  particularly 
effective  as  Galerkin  trial  functions  (or  basis  vectors).  Subsequently,  the  present  authors 
[2,5-7]  have  shown  that  the  terms  in  either  a  regular  or  a  singular  perturbation  expansion 
of  the  solution  are  also  effective  trial  functions.  In  particular,  we  have  demonstrated  that 
the  “reduced-basis”  solutions  can  be  useful  for  significantly  larger  values  of  the  expansion 
parameter  than  the  Taylor  series  or  singular  perturbation  solutions  on  which  they  are  based. 
A  treatment  of  the  reduced  basis  method  from  a  mathematical  point  of  view  is  given  by  Fink 
and  Rheinbolt  [3]. 

Some  general  observations  about  the  technique  are  the  following:  1)  In  many  perturbation 
problems,  much  effort  has  to  be  expended  to  compute  (analytically)  each  additional  term 
in  a  perturbation  expansion.  Through  the  use  of  the  proposed  hybrid  method,  the  known 
perturbation  terms  can  be  exploited  more  fully.  2)  Another  way  of  viewing  the  technique  is  to 
recognize  that  in  many  perturbation  expansions  the  functional  form  of  the  higher-order  terms 
can  be  well  approximated  by  a  linear  combination  of  the  lower-order  terms.  Thus,  much  of 
the  effect  of  the  higher  order  terms  may  be  included  by  applying  the  reduced  basis  technique 
to  a  small  number  of  lower  order  terms.  3)  Since  it  is  often  possible  to  develop  formal 
perturbation  expansions  for  the  solution  about  two  or  more  values  of  the  parameter,  e.g.  for 
small  or  large  values  of  the  parameter,  the  proposed  method  is  a  convenient  way  of  combining 
the  information  contained  in  these  different  expansions.  This  is  accomplished  by  allowing  the 
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set  of  basis  functions  to  include  terms  from  the  different  perturbation  expansions,  which  may 
be  either  regular  or  singular  expansions  (see  [7]  for  several  examples).  4)  While  the  use  of  a 
Taylor  series  expansion  is  frequently  limited  by  a  finite  radius  of  convergence,  the  proposed 
hybrid  method  can  sometimes  yield  good  results  even  well  outside  the  radius  of  convergence 
(see  [2]  as  well  as  our  first  example  here),  or  even  past  a  real  singularity  in  the  parameter  [6]. 
In  fact,  when  information  from  two  or  more  expansions  is  employed,  the  method  appears 
to  provide  meaningful  (and  often  very  accurate)  approximations  in  “intermediate”  regions 
of  parameter  values,  e.g.,  in  regions  where  the  parameter  is  neither  “small”  nor  “large”  (see 
especially  Geer  and  Andersen  [7j).  5)  The  same  general  technique  may  also  be  applied  in 
a  fully  numerical  sense  in  that  the  perturbation  functions  themselves  may  be  computed  in 
discretized  form.  6)  However  the  perturbation  functions  are  computed,  the  hybrid  solutions 
share  with  the  perturbation  solutions  the  property  that  the  accuracy  is  greatest  near  the 
expansion  points  and  diminishes  as  the  parameter  moves  away  from  these  points.  Usually 
a  good  indicator  of  accuracy  is  found  by  examining  the  difference  in  the  hybrid  solutions 
based  on  iV-1  terms  and  N  terms.  Thus,  it  is  possible  to  determine  the  range  of  parameter 
values  for  which  the  TV-term  expansion  is  valid. 

In  the  following  section,  we  describe  the  method  in  more  detail  and  then  apply  it  to  a 
simple  linear  boundary  value  problem  in  Section  3.  In  Section  4,  we  apply  the  method  to  the 
problem  of  determining  the  flow  of  a  slightly  compressible  fluid  past  a  circular  cylinder,  while 
in  Section  5  we  apply  it  to  the  problem  of  determining  the  shape  of  a  free  surface  induced 
by  the  presence  of  a  sink  above  the  surface.  Each  of  the  problems  in  sections  3-5  involves 
an  elliptic  boundary  value  problem,  yet  each  serves  as  a  model  problem  for  a  much  wider 
class  of  problems.  The  problem  of  Section  3  is  linear,  and  an  arbitrary  number  of  terms  in 
its  perturbation  solution  can  be  obtained  in  a  straightforward  manner.  Thus,  while  an  exact 
analytical  solution  is  not  known,  it  is  possible  to  perform  a  rather  detailed  analysis  of  both  the 
perturbation  and  hybrid  results,  including  some  explicit  expressions  for  some  of  the  hybrid 
solutions.  In  this  problem  the  radius  of  convergence  of  the  perturbation  solution  is  finite 
and  is  limited  by  singularities  which  lie  on  the  imaginary  axis  of  the  complex  parameter 
plane.  The  problem  of  Section  4  involves  a  nonlinear  partial  differential  equation,  whose 
perturbation  solution  is  tedious  to  obtain.  Its  radius  of  convergence  is  finite  and  appears  to 
be  limited  by  singularities  which  have  real  physical  significance.  The  perturbation  solution  of 
the  problem  of  Section  5  is  the  most  difficult  of  the  three  problems  to  compute  and  is  only  an 
asymptotic  (i.e. ,  a  nonconvergent)  expansion  [13].  For  each  of  these  problems  we  demonstrate 
that  the  hybrid  method  provides  a  much  more  useful  solution  than  the  perturbation  solution 
alone.  Some  general  comments  about  the  method  and  a  discussion  of  its  possible  application 
to  certain  problems  associated  with  domain  decomposition  are  discussed  in  Section  6. 
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2.  Description  of  the  Method 

The  method  we  wish  to  describe  is  a  two-step  hybrid  analysis  technique  based  on  the 
successive  use  of  perturbation  expansion  methods  and  the  classical  Bubnov- Galerkin  approx¬ 
imation  technique.  To  illustrate  the  general  ideas  of  the  method,  suppose  we  are  seeking  (an 
approximation  to)  the  solution  u  to  the  problem 

(2.1)  £(u,£)  =  0, 

where  £  is  some  partial  differential  operator  and  £  is  a  parameter.  (Although  we  restrict 
our  attention  to  two-dimensional  elliptic  boundary  value  problems  in  this  paper,  we  believe 
that  the  method  has  a  much  wider  range  of  applicability.  Hence,  we  formulate  the  method 
in  terms  of  a  general  operator  £.)  Here  (2.1)  holds  in  some  domain  V  ,  and,  in  addition, 
u  must  satisfy  certain  conditions  on  the  boundary  of  T>,  which  we  denote  by  dT>.  Without 
loss  of  generality,  we  can  assume  that  these  boundary  conditions  are  homogeneous  in  u. 

In  the  first  step  of  the  method,  we  generate  the  coordinate  functions  in  a  perturbation 
expansion  of  u  about  one  or  more  specific  values  of  the  parameter  e,  say  about  e  =  ep,  p  = 
1, 2, . . . ,  P.  In  the  second  step,  we  construct  new  approximate  solutions  consisting  of  sums  of 
some  of  these  perturbation  coordinate  functions,  each  multiplied  by  an  unknown  amplitude, 
and  then  determine  these  amplitudes  by  using  the  Bubnov-Galerkin  method. 

To  describe  this  idea  in  more  detail,  suppose  that  the  solution  to  (2.1)  can  be  expanded 
about  s  =  £p  into  a  series  of  the  form 

(2.2)  u  =  +  0(a^+1(£)), 

where  {^(e)}  is  an  appropriate  asymptotic  sequence  of  gauge  functions  and  each  u^'3'1  can  be 
determined  completely  by  standard  perturbation  methods  (see,  for  example,  Nayfeh  [8]).  By 
our  assumption  on  the  boundary  conditions  on  u,  each  u satisfies  homogeneous  boundary 
conditions  on  &D. 

A  subset  of  all  of  the  perturbation  functions  is  now  chosen  as  the  set  r .  coordinate 

functions  for  the  Bubnov-Galerkin  technique  and  approximations  u  for  u  arc  sought  in  the 
form 

N 

(2.3)  u=5>«fJiW(e), 

J=1 

where  the  (unknown)  parameters  {6>,f/(e)}  represent  the  amplitudes  of  the  coordinate  func¬ 
tions  iS3\  Here  each  u ^  is  one  of  the  perturbation  coordinate  functions  u^p'3\  and  hence 
u  satisfies  the  boundary  conditions  of  the  problem  for  any  choice  of  amplitudes  To 
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determine  these  amplitudes,  we  apply  the  Bubnov- Galerkin  technique  to  the  governing  equa¬ 
tion  (2.1).  Thus,  we  substitute  (2.3)  into  (2.1)  and  require  that  the  residual  be  orthogonal 
to  the  N  coordinate  functions  over  the  domain  T> ,  i.e., 


(2.4) 


1  <  k  <  N. 


Equations  (2.4)  represent  a  set  of  N  equations  for  the  N  unknown  amplitudes.  While  (2.4) 
must,  in  general,  be  solved  numerically,  solving  it  is  much  simpler  than  numerically  solving 

(2.1).  In  particular,  for  a  fixed  value  of  e,  the  solution  to  (2.4)  is  a  point  in  TV- dimensional 
space,  where  N  is  reasonably  small,  while  the  solution  of  (2.1)  is  a  continuous  function  of 
the  variables  x. 

We  should  note  that  this  particular  choice  of  coordinate  functions  overcomes  the  main 
drawback  of  the  Bubnov-Galerkin  method,  which  is  the  difficulty,  from  a  practical  point 
of  view,  of  selecting  a  small  number  of  good  coordinate  functions.  By  the  way  they  are 
constructed,  the  perturbation  coordinate  functions  are  (under  certain  assumptions)  elements 
of  a  set  of  functions  which  span  the  space  of  solutions  in  a  neighborhood  of  their  point  of 
generation.  Thus,  they  characterize  the  solution  u  in  that  neighborhood.  We  also  observe 
that,  in  many  applications,  the  functions  u^p'1^  are  determined  by  solving  a  set  of  linear 
equations,  even  though  the  original  operator  C  may  be  nonlinear. 

For  the  three  applications  discussed  in  this  paper  we  use  P  =  1  only. 


3.  A  Simple  Example 

To  illustrate  the  basic  features  of  the  hybrid  method,  we  consider  the  following  simple 
two-dimensional  example.  We  define  u(x,y,e)  as  the  solution  to  the  problem 

(3.1)  V2u  +  e  sin(x)  cos(y)  ux  +  2e  sin(x)  sin(y)  =  0,  ( x,y)eV , 


(3.2)  with  u  =  0  for  ( x,y)edV . 

Here  V  =  {(x,y)  :  0  <  x,y  <  7r}  and  dV  denotes  the  boundary  of  V.  In  (3.1),  the  sym¬ 
bol  V2  denotes  the  usual  two-dimensional  Laplacian  operator,  and  the  subscript  x  denotes 
differentiation  with  respect  to  x.  The  solution  seems  to  be  positive  over  the  whole  domain 
V  for  positive  values  of  e.  It  is  invariant  under  the  180°  rotation  x  — >  7r  —  x,  y  — >  7r  —  y. 
Equation  (3.1)  is  also  invariant  under  the  transformation  t  — >  —  e  and  u(x,y)  — ►  —u(x,n  —  y). 
Thus,  we  may  focus  our  attention  on  solutions  for  positive  e.  Figure  1  shows  surface  and 
contour  plots  of  the  solution  for  three  values  of  e.  For  large  values  of  e,  the  solution  exhibits 
a  number  of  boundary  layer  properties  (see  Section  6). 
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Step  One:  For  small  values  of  e,  we  construct  a  regular  perturbation  expansion  of  u  in 
the  form 

(3.3)  u  =  ]TV  u^3\x,y)  +  0(e^+1), 

3=1 

where  each  of  the  perturbation  coefficient  functions  u ^  is  independent  of  e.  To  determine 
these  functions,  we  substitute  (3.3)  into  (3.1)  and  (3.2)  and  equate  the  coefficients  of  like 
powers  of  e.  In  this  way,  we  are  led  to  the  following  system  of  equations,  from  which  the 
can  be  determined  recursively. 

(OA)  V72  (j) _  f  2  sin(x)  sin(y)  if  j  =  1 

'  (  —  sin(a:)  cos(y)  if  2  <  j  <  N, 

(3.5)  —  0  for  x  =  0,7r;  and  y  =  0, 7r  for  1  <  j  <  N. 

In  particular,  using  (3.4)  and  (3.5)  we  find 

=  sin(x)  sin(y),  u ^  =  (32)  sin(2z)  sin(2y), 

u(3)  =  (ik)  (1)  sin(3x)sin(3y)  +  Q)  sin(3x)sin(y) 

—  sin(x)  sin(3y)  —  sin(x)  sin(y)  , 

(3.6)  i  7 

I  sin(4x)  sin(4y) -f- (■———)  sinf4x)  sin(2w) 

V49152/  K  J  V76800/  K  J  K  y> 

—  (  sin(2x)  sin(4y)  —  ( — - — ^  sin(2x)  sin(2y). 

\  19200 /  v  J  K  \1920/  k  K  J 

It  is  easy  to  determine  the  general  form  of  each  u ^  and  to  construct  a  recurrence  relation 
for  the  numerical  coefficients  involved.  Thus,  many  more  terms  in  the  expansion  can  be 
computed  in  a  straightforward  manner. 

Step  Two:  In  the  second  step  of  the  hybrid  method,  we  use  the  perturbation  coordinate 
functions  { u &)}  obtained  in  Step  One  to  construct  new  approximate  solutions  u  in  the  form 

(3.7)  u=  £6jiN(e)u(j)(x,y), 

n=  1 

where  IV  is  a  relatively  small  integer  and  the  new  amplitudes  will  be  determined  using 

the  Galerkin  technique.  (We  note  that  the  u  defined  by  (3.7)  satisfy  the  boundary  condition 
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(3.2)  for  any  choice  of  the  because  of  the  conditions  (3.5)  on  the  functions  {u^}.) 

The  procedure  is  as  follows.  We  substitute  (3.7)  into  (3.1)  and  require  the  residual  to  be 
orthogonal  to  each  of  the  perturbation  coordinate  functions,  i.e., 


J  J  [^2/^  +  £  sin(x)  cos(y)ux  +  2e  sin(x)  sin(i/)]  u W  dx  dy  —  0, 


which  may  be  written  in  the  form 


(3.8) 


N 


T  £  dk,j]  —  £  bj, 

j'= i 


1  <  k  <  N, 


where 


Ck,j  =  J  J  [V2u^]  dx  dy,  dkj  =  J  J  [sin(x)  cos(y)  u^]  dx  dy, 


(3.9) 


bj  =  —  J  J  2sin(x)sin (y)u^dxdy. 


For  this  simple  example,  since  the  original  problem  is  linear,  the  equations  (3.8)  to  determine 
the  are  linear.  Further,  since  the  integrals  in  (3.9)  can  be  evaluated  exactly,  it  is 

possible  in  this  case  to  evaluate  the  amplitudes  as  rational  functions  of  e.  However,  we  remark 
that  for  most  applications  (especially  nonlinear  ones)  the  evaluation  of  the  amplitudes  must 
be  carried  out  numerically,  typically  by  an  iterative  process  for  gradually  increasing  values 
of  the  parameter  c. 

For  N  =  2,  we  use  the  expressions  (3.6)  to  evaluate  the  coefficients  {ck,j},  {dk,j}  and 
explicitly  and  find 


(3.10) 


<$1,2  = 


1  +  c2/128’ 


32,2 


1  +  c2/128‘ 


In  a  similar  manner,  using  the  symbolic  manipulation  system  Mathematica  [15],  for  N  =  3 
we  find 

£2 

(3.11)  «,.3  =  e,  «m=  i+£2/60.  *3.3  =  ^=,,, 

while  for  N  =  4 


56524800  c  +  1216128  c3 
56524800  +  1216128  c2  +  2141c4’  2,4  “  £  M’ 


(3.12) 


56524800  c3 

56524800  f  1216128  c2  +  2141  c4  ’  4,4  “  £  3,4 
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Using  Mathematica,  we  have  obtained  explicit  expressions  for  the  6j^  for  N  <  8.  In  each 
case,  we  find  that 

(3.13)  8hN  =  ej  +  0(eN+1)  as  e  ->  0,  for  1  <  j  <  N. 

Hence  we  have  verified  that  our  N-term  hybrid  solution  (3.7)  agrees  with  the  iV-term  per¬ 
turbation  solution  (3.3)  as  e  — *  0  through  N  =  8. 

In  Figure  2,  we  have  plotted  the  relative  Z^-error  between  numerical  solutions  to  the 
problem  (3.1 )— (3.2)  and  either  hybrid  solutions  (3.7)  or  perturbation  solutions  (3.3)  for  dif¬ 
ferent  values  of  N.  Figure  2  demonstrates  that:  (1)  there  is  a  finite  radius  of  convergence 
for  the  perturbation  solutions;  (2)  the  hybrid  solutions  are  more  accurate  than  the  corre¬ 
sponding  perturbation  solutions,  especially  for  larger  values  of  e;  and  (3)  for  a  fixed  N  the 
accuracy  of  the  hybrid  solutions  gradually  deteriorates  as  £  increases.  We  discuss  this  exam¬ 
ple  further  in  Section  6.  The  numerical  solutions  were  obtained  using  a  simple  second  order 
finite  difference  method  (coupled  with  SOR)  on  a  41  by  41  uniform  grid.  For  large  values 
of  £,  the  numerical  solutions  are  more  accurate  than  the  perturbation  and  hybrid  solutions, 
although  for  small  e  the  reverse  is  true. 


4.  Circular  Cylinder  in  Slightly  Compressible  Flow 


We  now  apply  our  hybrid  method  to  a  problem  involving  a  nonlinear  partial  differential 
equation.  As  a  model  of  this  type  of  problem,  we  consider  the  problem  of  determining  the 
steady,  two-dimensional  flow  of  an  inviscid,  compressible,  perfect  gas  past  a  circular  cylinder 
of  unit  radius  without  circulation  (see,  e.g.,  Van  Dyke  [10]).  If  we  introduce  the  usual  polar 
coordinates  (r,  6)  and  express  the  fluid  velocity  q  in  terms  of  a  potential  function  ip  as 


U  V{  (r  +  r  cos(0)  +  v?(r,0,  M)| , 


we  find  that,  for  r  >  1,  <p  satisfies 


=  VV-(|)M2{[(l-r-2)cos(0)  +  ¥>r]Qr 


(4-1) 

+  [— (r  1  +  r  3)  sin($)  +  r  2  ipg]Qg  +  -  l)  Q  V2<p}  =  0, 

where 

Q  =  r-4  —  2  r~2  cos(20)  +  2  cos(0)  (1  —  r-2)  <pr  —  2  sin($)  (r_1  +  r~3)  ipg 

+  ¥>2  +  r-2¥>2, 


with  Ifir  =  0  on  r  =  1,  and  <p  —  0(r  ')  as  r  — >  oo.  Here  M  =  U/c  is  the  free  stream  Mach 
number,  c  is  the  speed  of  sound  in  the  gas,  and  7  is  the  adiabatic  ratio.  For  small  values 
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of  M  the  flow  streamlines  are  shown  in  Figure  3.  For  subsonic  flow,  the  maximum  velocity 
occurs  at  the  surface  of  the  cylinder  at  9  —  7r/2  and  9  —  3tt/2.  We  now  apply  the  hybrid 
method  as  outlined  in  Section  2  to  this  (nonlinear)  problem,  where  M  plays  the  role  of  the 
parameter  z. 

Step  One:  For  small  values  of  M,  we  use  the  regular  perturbation  method  to  obtain 

N 

(4.2)  <p(r,  9,  M)  —  M2j  <p(J\r,  9)  +  0(M2N+2)  as  M  ->  0, 

j=i 

where  each  satisfies  a  Poisson  equation  in  the  region  exterior  to  the  cylinder,  satisfies  the 
condition  p[^  =  0  on  r  =  1,  and  is  0(r_1)  as  r  — >  oo.  The  general  form  of  the  perturbation 
coordinate  functions  {p^}  is  given  by 

(4.3)  <P(;)  =  H /j.fcO")  cos[(2fc  + 1)0],  j  >  1. 

k= o 


In  particular,  we  find: 


Van  Dyke  and  Guttmann  [12]  have  computed  29  terms  in  the  expansion  (4.2)  using  a  FOR¬ 
TRAN  program  and  have  presented  a  rather  detailed  analysis  of  the  convergence  of  the 
series  as  IV  ->  oo.  We  have  used  the  symbolic  computation  system  Mathematica  [15]  to 
compute  5  of  these  terms  in  exact  rational  arithmetic,  carrying  7  as  a  parameter,  and  11 
terms  for  the  special  case  when  7  =  7/5.  Although  the  general  form  of  each  p^  is  known 
and  only  certain  numerical  coefficients  need  to  be  determined,  the  amount  of  computational 
“labor”  (using  either  purely  numerical  or  symbolic  computation)  increases  significantly  as 
N  increases.  Consequently,  it  is  desirable  to  obtain  as  much  information  as  possible  about 
tl  3  solution  from  the  first  few  perturbation  coefficient  functions  {p^}.  This  we  do  in  Step 
Two  below. 

Step  Two:  We  now  seek  new  approximate  solutions  (p  in  the  form 
(4.5)  if{r,e,M)  =  YI6hN(M)  VW(r,«), 

J=  1 
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where  the  new  “amplitudes”  {£j,A/-}  must  be  determined.  To  determine  the  {^jv},  we  sub¬ 
stitute  (4.5)  into  the  governing  differential  equation  (4.1)  and  require  that  the  residual  be 
orthogonal  to  each  of  the  perturbation  coordinate  functions  ip^k\  i . e. , 

/■2ir  r  oo 

(4.6)  J  j  £(<p,M)<pWrdrd&  =  0,  1  <  k  <  N. 

Equations  (4.6)  are  a  system  of  N  (cubically  nonlinear)  equations  to  determine  the  N  am¬ 
plitudes  {6j'n(M)}.  They  must,  in  general,  be  solved  numerically.  However,  this  is  straight¬ 
forward  to  do  using  a  standard  method  such  as  Newton’s  method.  In  particular,  by  starting 
at  “small”  values  of  M,  where  we  expect  8hw  m  M2\  and  then  proceeding  to  larger  values 
of  M,  the  solution  of  (4.6)  can  be  obtained  in  an  efficient  manner. 

For  the  special  case  when  N  =  1,  we  let  (p  —  £lfl  where  c,^1)  is  given  by  (4.3)  and 
(4.4)  with  j  =  1,  and  we  find  that  equation  (4.6)  yields  the  following  cubic  equation  for 

8  =  5U: 

(4.7)  8-  M2[l  +  c18  +  c2  62  +  c383}  =  0, 
where 

cx  =  (4736  +  1111 7)/7532,  c2  =  (147139  +  158248  ^)/1084608, 

(4.8) 

c3  =  (9182  +  815931  t)/23861376. 

From  (4.7)  and  (4.8),  we  see  that 

(4.9)  8  =  M2  +  ci  M4  +  0(M6)  asM->0, 

and  hence  our  hybrid  solution  reproduces  the  first  term  of  the  perturbation  solution  as 
M  —>  0.  Also,  the  cubic  equation  (4.7)  has  one  negative  real  root  and  two  positive  real  roots 
for  0  <  M  <  Mc  %  0.5389247.  For  M  >  Mc,  there  are  no  positive  real  roots  of  this  equation. 
In  Section  6  we  discuss  a  possible  interpretation  of  Mc  in  relation  to  the  convergence  of  the 
series  (4.2). 

In  a  similar  manner,  for  N  =  2, 3, 4,  and  5,  we  have  used  Mathematica  to  perform  the 
integrations  appearing  in  (4.6)  and  have  expressed  the  resulting  equations  in  exact  rational 
arithmetic.  We  then  expressed  the  solutions  as  Taylor  series  in  M 2  and,  in  each  case, 

we  found  that 

(4.10)  6jtN  =  M2’  +  0(M2N+2)  as  M  ->  0,  for  1  <  j  <  N  . 

Thus,  for  these  values  of  N  we  have  verified  that  our  hybrid  solution  (4.5)  agrees  with  the 
first  N  terms  in  the  perturbation  expansion  (4.2)  as  M  — >  0.  In  addition,  for  each  of  these 
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values  of  N,  we  have  also  estimated  the  quantity  Mc,  i.e.,  the  value  of  M  for  which,  if 
M  >  Mc,  me  system  (4.6)  does  not  appear  to  have  physically  meaningful  solutions.  These 
values  of  Mc  are  summarized  in  Table  1. 

As  an  application  of  the  results  we  have  obtained,  we  estimate  the  critical  Mach  number 
M,  associated  with  this  flow.  Here  M,  is  defined  as  the  value  of  the  free-stream  Mach  number 
at  which  the  flow  locally  becomes  sonic  for  the  first  time  and  is  determined  as  the  solution 
of  the  equation 

Ml)  M.j=2/[(7  +  1)9L,-7  +  1). 

Here  qmax  —  2  —  dp>/39  evaluated  at  r  =  1,  and  6  =  7r/2  is  the  maximum  surface  speed. 
Van  Dyke  and  Guttmann  [12]  have  used  29  terms  in  the  expansion  (4.2),  nlong  with  a 
battery  of  numerical  techniques  (see  also  Andersen  and  Geer  [1]  and  Van  Dyke  [11]),  to 
obtain  the  estimate  M,  =  0.39823780  ±  0.0000001.  In  contrast,  using  our  one-term  hybrid 
approximation  (p,  we  find  that  our  approximation  qmax  for  qmax  is  given  by  qm&x  =  2  + 
(7/6)  #1,1.  Using  equations  (4.7)  and  (4.8)  to  compute  for  a  particular  value  of  M,  we 
find  from  (4.11)  that,  for  7  =  7/5,  our  one-term  hybrid  solution  yields  the  approximation 
M,  =  0.407257966.  This  value  differs  from  Van  Dyke  and  Guttmann’s  estimate  by  only 
about  2.3%. 

In  a  similar  manner,  we  find  approximations  to  both  Mc  and  M.  using  N  —  2,3,4,  and  5 
terms  in  our  hybrid  approximation  (4.5).  The  estimates  we  obtain,  along  with  those  obtained 
by  using  N  terms  in  tne  perturbation  expansion  (4.2),  are  summarized  in  Table  1. 

TABLE  1 

Approximations  to  M, 


N 

Perturbation 

Hybrid 

Me 

1 

0.420943 

0.407258 

0.5389 

2 

0.409239 

0.401704 

0.4990 

3 

0.404577 

0.399699 

0.4740 

4 

0.402274 

0.398948 

0.4606 

5 

0.400979 

0.4513 

6 

0.400187 

7 

0.399671 

8 

0.399319 

9 

0.399071 

10 

0.398891 

11 

0.398757 

10 


In  Section  6,  we  discuss  certain  implications  of  the  values  which  appear  in  Table  1.  For 
the  present,  we  just  note  that  the  hybrid  approximation  to  M*  based  on  N  terms  has  about 
the  same  accuracy  as  the  perturbation  approximation  based  on  2N  terms.  We  have  not  tried 
to  push  our  technique  nearly  far  enough  to  begin  to  compete  with  Van  Dyke’s  results.  We 
only  observe  that  with  a  small  number  of  terms  we  have  a  better  approximation. 

5.  Sink  or  Source  Near  a  Free  Surface 

We  now  wish  to  apply  our  hybrid  method  to  a  problem  for  which  the  nonlinearity  enters 
through  the  boundary  conditions;  in  particular,  a  problem  in  which  the  boundary  itself  is 
unknown  and  must  be  determined  as  part  of  the  solution  to  the  problem.  As  a  model  of 
this  type  of  problem,  we  consider  the  problem  of  determining  the  steady,  two-dimensional 
shape  of  a  free  surface  when  a  sink  or  source  is  placed  above  the  surface.  This  problem  has 
been  studied  using  a  variety  of  numerical  techniques  by  several  investigators,  e.g.,  Tuck  and 
Vanden-Broeck  [9]  and  Vanden-Broeck  and  Keller  [14],  and  also  from  a  perturbation  point 
of  view  by  Vanden-Broeck,  Schwartz  and  Tuck  [13]. 

We  consider  a  sink  (or  source)  of  strength  Q  located  a  distance  h  above  the  undisturbed 
height  of  a  free  surface.  We  introduce  a  rectangular  coordinate  system  ( h  x,  hy ),  with  gravity 
acting  in  the  negative  y-direction  and  with  the  origin  a  distance  h  below  the  sink  (see  Figure 
4).  We  let  the  velocity  potential  be  $  =  Q  {(4")  log[x2  -f  (y  —  l)2]  +  <^(2,  y)}  and  denote  the 
elevation  of  the  free  surface  by  y  =  77(2).  Then  the  problem  to  determine  and  77  becomes 

(5.1)  VV  =  0,  y  >  77(2), 
with 

(5.2)  =  (A)  =  0 

and 

1  2 

(5.3)  B3(<p,y)  =  77  -  £077(1  +  t)2)~s  —  g  (1  +  7 )2)  ipx  +  x2  +  ^_1)2  =0 

on  y  =  t](x).  Here  =  d/dx ,  e  =  Q2/2g/i3,  and  a  =  2'yh/pQ2,  where  g  is  the  acceleration 
due  to  gravity,  p  is  the  density  of  the  fluid,  and  7  is  the  surface  tension  coefficient.  Equation 
(5.2)  is  just  the  kinematic  boundary  condition  on  the  free  surface,  while  equation  (5.3)  follows 
from  Bernoulli’s  equation  and  the  condition  on  the  jump  in  pressure  at  the  free  surface  due  to 
surface  tension.  We  are  interested  in  finding  approximations  to  the  solutions  ip  =  ip(x,y,e ) 
and  77  —  77(2,  e)  to  equations  (5. 1 )— (5.3)  for  small  values  of  e. 
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Step  One:  For  small  values  of  e,  we  look  for  regular  perturbation  solutions  tor  if  and  r/ 
in  the  form 

N- 1  N 

(5.4)  ip  =  ^2  e3  <pW  +  0(eN),  r)  =  T]j  +  0(eN+- )  as  e  — »  0. 

3=0  j=l 

Using  the  expressions  (5.4)  in  equations  (5.1 )— (5 .3) ,  we  find  that  the  terms  which  are  0(1) 
in  e  in  equations  (5.1)  and  (5.2)  yield  the  result 

(5.5)  if{0)  =  log(x2  +  (y  +  l)2). 

Then  the  terms  which  are  0(e)  in  equation  (5.3)  yield  the  expression 

(5.6)  7]i  —  x2  (1  +  a:2)-2. 

In  a  similar  manner,  the  terms  which  are  0(e)  in  equations  (5.1)  and  (5.2)  yield  the  result 
7)  (i)  _  J_  [3(z2  -  (y  +  l)2)  2 (y  +  l)\(y  +  l)2  -  3x2]l 


87T3  [(x2  +  (y  +  l)2)2  (x2  +  (y  +  l)2)3 

and  the  terms  which  are  0(e2)  in  (5.3)  yield  the  expression 


3 

x2(l  —  6  x2  +  x4) 

2  a 
_ 

T  -  8x2  +  3x4' 

2tt4 

(1  +  x2)5  j 

'  9 

(1  +  x2)4 

Continuing  this  procedure  with  the  aid  of  Mathematica,  we  find  that  the  general  form  of 
the  perturbation  coefficient  functions  ip ^  and  Tjj ,  j  >  1,  is 

(5.9)  ip{3\x,y)  =  (l/^+l)J2cJlk^k\x,y), 


Vj{x)  =  Pi.fcC^2)  C1  +  x2)  (3J  1  k)- 

k=0 

Here  each  ch is  a  constant,  each  py*,  is  a  polynomial  in  x2,  and  the  functions  {tp^}  are 
defined  by 

V>(0)  =  log[x2  +  (i/  +  l)2]5,  =  (d/dy)^k-V  for  Jb  >  1. 

In  particular,  we  find 

.  .  x2  (45  —  1909  x2  +  6890  x4  —  2786  x6  +  9  x8  —  9  x10) 

”3(l)  =  - 32  7T6  ( 1  +  x2)® - 

a  (12  -  853  X2  +  4008  x4  -  2658  x6  +  148  x8  -  x10) 

4  7T4  (  1  +  X2)7 

24  a2  (2  —  33  x2  +  40  x4  -  5x6) 

7T2(1+X2)6 
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(27  -  807r2a)/384, 


C2,2  =  -c2, 3  =  (9  -  87r2a)/128,  c2,4  = 

c2,5  =  (3  -  47T2a)/64,  c2,6  =  3/640. 

From  the  analysis  of  this  perturbation  series  with  zero  surface  tension  by  Vanden-Broeck 
et  al  [13]  (see  Section  6  for  more  discussion  of  their  results),  our  formal  series  expansion  of 
the  solution  for  this  example  does  not  converge  for  any  non-zero  value  of  e.  However,  despite 
this  lack  of  convergence,  we  proceed  to  Step  Two  of  our  hybrid  method. 

Step  Two:  Using  the  expressions  above  for  the  perturbation  coordinate  functions  ip^\x,y) 
and  T]j(x),  we  look  for  new  approximate  solutions  in  the  form: 

N  N 

(5.10)  v(x,e)  =  X^nOO^z),  ¥>(x,y,e)  = 

j=i  j=i 

where  the  new  “amplitudes”  must  be  determined.  (We  note  that  <p  satisfies  V2<p  =  0 

for  any  choice  of  the  amplitudes  {<5j,w}.)  To  determine  these  amplitudes,  we  substitute  the 
expressions  (5.10)  into  the  boundary  conditions  (5.2)  and  (5.3)  and  require  that  the  residuals 
be  orthogonal  to  appropriate  sets  of  test  functions  {a*}  and  {/3k},  he., 

/oo  roo 

Bx{(p,f})ak{x)dx  —  0,  /  fik(x)  dx  —  0,  1  <  k  <  N. 

-oo  J- oo 

Once  the  test  functions  have  been  selected,  equations  (5.11)  are  a  set  of  2N  equations 
for  the  2  N  unknowns  {6j,n,  j  =  1,2,  ...,2  N}.  For  simplicity  we  select  as  test  functions 
<*k  =  Pk  —  Vk(x)-  The  resulting  hybrid  solutions  are  discussed  in  the  next  section. 

6.  Discussion 

The  hybrid  perturbation-Galerkin  method,  as  we  have  described  it  here,  is  a  semi- 
analytical,  semi- numerical  technique  which  appears  to  have  the  potential  of  being  a  useful 
tool  both  to  complement  and  to  supplement  existing  standard  methods  of  analysis.  It  is 
semi-analytical  in  the  sense  that  some  of  the  analytical  structure  of  the  solution  is  deter¬ 
mined  by  first  constructing  one  or  more  perturbation  approximations  to  the  solution.  It  is 
semi-numerical  in  that  new  amplitudes  of  the  perturbation  coordinate  functions  typically  are 
determined  numerically  by  solving  standard  Galerkin  equations  associated  with  the  prob¬ 
lem.  The  method  seems  to  have  considerable  potential  to  be  an  effective  problem-solving 
tool  because  the  perturbation  coordinate  functions  appear  to  be  very  effective  trial  functions 
in  the  Galerkin  approximation.  Intuitively  this  is  because,  by  the  way  they  are  constructed, 
they  describe  the  basic  nature  of  the  solution,  at  least  for  parameter  values  near  their  point 
of  generation. 
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In  our  first  example  (Section  3),  we  analyzed  a  simple,  linear  boundary-value  problem 
on  a  bounded  domain,  for  which  an  arbitrary  number  of  terms  in  the  regular  perturbation 
series  can  be  computed  in  a  straightforward  manner.  By  examining  the  singularities  of  the 
amplitudes  {<5j,/y}  in  our  hybrid  solution,  it  appears  that  the  convergence  of  the  perturbation 
series  is  limited  by  a  pair  of  complex  conjugate  singularities  located  near  or  at  ±6.70±  (This 
follows  by  analogy  from  the  analysis  of  several  similar  problems,  involving  ordinary  differ¬ 
ential  equations,  for  which  the  singularities  of  the  {£_;,//}  in  the  complex  £-plane  converged 
to  the  singularities  of  the  actual  solution  as  the  number  of  terms  in  the  approximation  in¬ 
creased.  See  especially  [2],  In  fact  the  location  of  the  poles  of  the  which  are  closest 

to  the  origin  form  the  sequence  ±7.75  i,  ±7.15  i,  ±6.803  i,  ±6.721  i,  ±6.7048  i,  ±6.70283  i, 
Another  pair  of  singularities  is  indicated  near  e  ~  ±9t  .)  These  singularities  appear  to  have 
no  immediate  interpretation,  but  nonetheless  limit  the  direct  usefulness  of  the  perturbation 
expansion. 

In  Figure  2,  we  have  plotted  the  relative  T2-error  between  numerical  solutions  to  problem 
(3.1 )— (3.2)  and  perturbation  or  hybrid  approximations.  The  figure  clearly  illustrates  that 
the  hybrid  method  allows  us  to  obtain  useful  information  about  the  solution  well  beyond 
the  radius  of  convergence  of  the  perturbation  solution.  Also,  there  is  no  indication  that  the 
range  of  convergence  of  the  hybrid  solution  is  limited. 

In  Figure  6,  we  have  plotted  the  values  of  the  solution  at  the  center  value  x  —  y  —  7r/2  as 
calculated  numerically  and  by  the  use  of  the  hybrid  technique,  where  the  number  of  terms 
( N )  ranges  from  1  to  8.  It  is  interesting  to  note  that  our  hybrid  solutions  appear  to  form  an 
alternating  sequence  which  brackets  the  true  center  point  solution.  Figure  1  shows  that  the 
solutions  have  their  peak  “elevations”  at  the  center  point. 

The  contour  plots  shown  in  Figure  7  demonstrate  that  higher-order  perturbation  func¬ 
tions  can  be  very  similar  to  one  another.  This  suggests  that  the  perturbation  functions  which 
are  of  higher  order  than  these  shown  may  be  fairly  well  approximated  by  linear  combinations 
of  the  functions  shown  (see  observation  (2)  in  the  introduction).  Figure  7  also  demonstrates 
that  the  perturbation  functions  may  have  higher  symmetries  (group  order  =  4)  than  the  so¬ 
lutions  (group  order  =  2)  which  they  are  used  to  approximate.  Such  additional  symmetries 
may  be  exploited  in  the  steps  of  the  hybrid  method  which  are  computed  numerically.  The 
symmetries  of  the  perturbation  functions  follow  a  regular  pattern  which  can  be  determined 
a  priori. 

From  a  singular  perturbation  point  of  view,  as  e  becomes  large  the  solution  to  problem 
(3.1 )— (3.2)  takes  on  a  somewhat  complicated  form.  The  leading  term  in  the  outer  expansion 
of  u  is  2(7r  -  x)tan(y)  for  0  <  x  <  n,  0  <  y  <  7t/2,  and  is  -2xtan(y)  for  0  <  x  < 
7r,  7r/2  <  y  <  7T.  Boundary  layers  develop  along  x  =  0  for  0  <  y  <  7r/2  and  along  x  =  7r 
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for  7r/2  <  y  <  tt.  In  addition,  an  internal  layer  forms  along  y  —  7r/2,  as  well  as  additional 
boundary  layers  at  (0,7r/2)  and  (m ,  7r /2) .  From  the  contour  plots  shown  in  Figure  1,  it  may 
be  seen  that  the  boundary  and  internal  layers  are  beginning  to  be  formed  by  the  hybrid 
approximations  as  e  increases,  even  though  these  approximations  are  based  on  the  small  e 
expansion  of  the  solution,  where  no  such  layers  are  evident. 

This  last  observation  suggests  that  the  hybrid  method  might  be  a  useful  tool  for  the 
general  problem  of  determining  an  appropriate  decomposition  of  the  domain  for  a  purely 
numerical  solution  to  a  specific  problem.  In  particular,  suppose  that  some  “singular”  behav¬ 
ior  of  the  solution  is  anticipated  for  large  values  of  a  parameter  appearing  in  the  problem 
formulation.  To  help  determine  the  location  and,  to  some  extent,  the  nature  of  this  singular 
behavior,  the  hybrid  method  could  be  applied  in  the  following  way.  First,  a  (regular)  pertur¬ 
bation  expansion  of  the  solution  is  constructed  for  small  values  of  the  parameter.  The  hybrid 
method  is  then  applied,  using  the  small  parameter  perturbation  coordinate  functions,  and 
then  the  behavior  of  the  hybrid  approximation  is  noted  as  the  parameter  value  is  increased. 
From  the  example  of  Section  3,  as  well  as  several  other  examples  we  have  examined  in  some 
detail,  it  appears  that  the  hybrid  solution  simulates  at  least  some  of  the  singular  behavior 
of  the  exact  solution  as  the  parameter  value  is  increased,  and  that  this  simulation  becomes 
more  accurate  as  the  number  of  terms  in  the  approximation  is  increased. 

In  our  second  example  (Section  4),  which  involves  an  infinite  domain,  the  regular  per¬ 
turbation  expansion  of  ip  is  straightforward,  although  tedious,  to  compute.  Frankl  and 
Keldysh  [4]  have  proved  that  y?  is  analytic  in  M 2  about  M2  =  0  and  hence  the  pertur¬ 
bation  expansion  (4.2)  is  actually  the  Taylor  series  expansion  of  the  exact  solution.  They 
did  not,  however,  determine  the  radius  of  convergence  Mc  of  the  series.  Through  a  care¬ 
ful  analysis  of  the  first  24  terms  in  the  series,  Van  Dyke  and  Guttmann  [12]  report  that 
Mc  =  0.402667605  ±  0.000000005,  which  is  about  1.11%  greater  than  their  estimate  of  Mt. 
If  this  were  indeed  the  case,  it  would  imply  that  shock  free  solutions  exist  for  a  continu¬ 
ous  range  of  values  of  M  above  the  critical  Mach  number.  However,  Van  Dyke  [private 
communication]  now  believes  this  result  to  be  in  error  and  that  Mc  =  Mt. 

Although  Van  Dyke  and  Guttmann  were  unable  to  determine  the  nature  of  the  singularity 
which  limits  the  convergence  of  the  series,  we  can  safely  assume  that  the  singularity  is  related 
to  the  formation  of  a  shock  and  hence  represents  a  real  singularity  with  a  definite  physical 
interpretation.  The  numbers  Mc  presented  in  Table  I  are  estimates  for  the  value  of  M  above 
which  no  physically  meaningful  hybrid  solution  to  the  problem  exists.  In  this  sense,  they  can 
be  thought  of  as  estimates  of  the  radius  of  convergence  of  the  series  (4.2).  If  more  of  these 
estimates  were  computed  and  if  it  could  be  shown  that  they  converge  to  the  same  limit  as 
the  estimates  A?*,  this  would  lend  support  to  the  conjecture  that  Mc  =  M,. 


The  solution  to  our  third  example  (Section  5)  involves  two  unknown  functions  (i.e. , 
the  potential  function  and  the  shape  of  the  free  surface),  and  the  associated  perturbation 
expansions  are  somewhat  more  difficult  to  compute  than  the  expansions  associated  with 
the  first  two  examples.  However,  the  general  form  of  the  terms  in  each  expansion  can  be 
determined  and,  in  principle,  only  certain  numerical  constants  need  to  be  computed.  Using 
a  complex  variable  formulation  of  this  problem  (with  zero  surface  tension)  Vanden-Broeck  et 
al  [13]  numerically  computed  these  constants  in  the  first  34  terms  in  each  expansion.  They 
analyzed  the  resulting  series  for  the  free  surface  and  found  it  to  be  of  exponential-integral 
character  (i.e.,  the  coefficient  of  £n  behaves  like  n!  for  large  n).  Consequently,  the  series 
diverges  for  all  nonzero  values  of  e.  They  use  a  number  of  techniques  to  “sum”  this  divergent 
series  and  show  that  this  analysis  leads  to  a  free  surface  shape  with  a  jump  discontinuity 
due  to  the  location  of  a  branch  cut.  One  such  free  surface  profile  corresponding  to  our 
£  =  7r2/5  =  1.974  is  indicated  by  a  dashed  line  in  Figure  5.  Our  hybrid  solution  with  N  =  2 
is  shown  by  a  solid  line.  Vanden-Broeck  et  al  indicate  that  the  jump  discontinuity  in  their 
solution  could  be  removed  a  applying  a  certain  iterative  procedure,  but  the  corresponding 
profile  was  not  displayed.  Using  a  numerical,  series  truncation  method  on  another  complex 
variable  formulation  of  the  problem,  Tuck  and  Vanden-Broeck  [9]  found  a  cusp  like  solution, 
corresponding  in  our  notation  to  a  value  of  £  =  £c  %  6.311.  Vanden-Broeck  and  Keller 
[14],  using  a  similar  method,  treated  the  same  problem  we  are  considering,  except  that  a 
horizontal  fixed  surface  is  present  at  a  finite  distance  above  the  sink.  They  found  solutions 
corresponding  to  our  problem  for  values  of  £  larger  than  ec,  but  did  not  find  any  steady 
solutions  for  0  <  £  <  £c  .  They  indicate,  however,  that  solutions  with  waves  may  exist  for  £ 
in  the  range  0  <  e  <  £c. 

The  hybrid  approximations  presented  in  Section  5  appear  to  give  “physically  realistic” 
approximations  to  a  steady  solution  to  the  problem,  as  indicated  in  Figures  4  and  5.  In 
[2]  we  applied  our  hybrid  method  to  an  eigenvalue  problem  associated  with  an  anharmonic 
oscillator.  For  this  problem,  the  perturbation  series  also  has  a  zero  radius  of  convergence, 
but  the  hybrid  approximations  converge  monotonically  to  the  true  solution  as  N  increases. 
Based  upon  this  experience,  it  is  certainly  possible  that  our  hybrid  approximations  to  the 
free  surface  for  values  of  £  in  the  range  0  <  e  <  tc  are  converging  to  a  steady  state  solution. 
Obviously,  this  difficult  problem  needs  much  more  investigation,  e.g.,  either  the  development 
of  a  numerical  method  which  will  yield  the  steady  solution  for  £  in  the  range  0  <  £  <  £c,  or 
a  proof  that  no  such  solution  exists  in  this  range  of  parameter  values. 
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FIGURE  CAPTIONS 


Figure  1.  Surface  and  contour  plots  for  8-term  hybrid  and  numerical  solutions  of  the  simple 
partial  differential  equation  of  Equations  (3.1 )— (3 .2) .  The  hybrid  solutions  shown  for  e  =  1 
and  e  =  5  are  highly  accurate.  Some  irregularities  may  be  seen  in  the  relatively  flat  portions 
of  the  hybrid  solution  for  e  =  20  which  do  not  appear  in  the  numerical  solution  for  e  =  20. 
Nevertheless,  the  relative  Z^-error  is  less  than  1%. 

Figure  2.  Relative  Z-2-errors  of  perturbation  and  hybrid  solutions  of  the  simple  PDE  problem 
as  a  function  of  the  parameter  e.  The  number  of  terms  N  ranges  from  1  to  8.  The  log-log 
plots  of  the  perturbation  errors  fall  nearly  on  straight  lines  which  intersect  one  another  near 
the  radius  of  convergence,  R  —  6.7026.  The  hybrid  errors  for  given  N  and  e  are  significantly 
less  than  the  perturbation  errors  for  the  same  N  and  e  and  show  no  signs  of  divergence. 

Figure  3.  Streamlines  for  two-dimensional  flow  around  a  cylinder  for  Mach  number  M  =  0, 
the  low  velocity  limit  of  the  solutions  to  Equation  (4.1). 

Figure  4.  Free  surface  profile  for  the  two-dimensional  flow  induced  by  a  sink  (black  dot) 
above  a  liquid  surface  (see  Equations  (5.1)-(5.3)).  The  effects  of  surface  tension  are  ignored. 
The  horizontal  and  vertical  scales  differ. 

Figure  5.  Free  surface  profiles  for  e  =  7r2/5  with  zero  surface  tension  as  computed  by  the 
hybrid  method  (solid  line)  and  as  computed  by  Vanden-Broeck  et  al  (dashed  line). 

Figure  6.  Maximum  values  of  the  solution  (x  =  y  =  7r/2)  of  the  simple  PDE  problem 
of  Section  3  as  computed  numerically  (dots)  and  by  the  hybrid  method  (solid  lines)  using 
from  N  =  1  to  N  =  8  terms.  The  results  for  N  =  8  are  useful  well  beyond  the  radius  of 
convergence,  R,  of  the  series  expansion. 

Figure  7.  Contour  plots  of  the  first  eight  perturbation  functions  for  the  simple  PDE  problem 
of  Section  3  demonstrate  that  the  higher  order  perturbation  functions  can  be  very  similar  to 
one  another  and  can  have  higher  symmetries  than  the  solutions  they  are  used  to  approximate. 
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